home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / GNU / gnuplot.lha / gnuplot / src / fit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-08-03  |  34.4 KB  |  1,249 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: fit.c 1.3 1993/08/03 07:14:17 cg Exp cg $";
  3. #endif
  4.  
  5. /*
  6.  *    Nonlinear least squares fit according to the
  7.  *    Marquardt-Levenberg-algorithm
  8.  *
  9.  *    added as Patch to Gnuplot (v3.2 and higher)
  10.  *    by Carsten Grammes
  11.  *    Experimental Physics, University of Saarbruecken, Germany
  12.  *
  13.  *    Internet address: cagr@rz.uni-sb.de
  14.  *
  15.  *    Copyright of this module:   Carsten Grammes, 1993
  16.  *
  17.  *    Permission to use, copy, and distribute this software and its
  18.  *    documentation for any purpose with or without fee is hereby granted,
  19.  *    provided that the above copyright notice appear in all copies and
  20.  *    that both that copyright notice and this permission notice appear
  21.  *    in supporting documentation.
  22.  *
  23.  *    This software is provided "as is" without express or implied warranty.
  24.  *
  25.  *    930726:     Recoding of the Unix-like raw console I/O routines by:
  26.  *            Michele Marziani (marziani@ferrara.infn.it)
  27.  */
  28.  
  29.  
  30. #define FIT_MAIN
  31.  
  32. #include <stdlib.h>
  33. #include <stdio.h>
  34.  
  35. #include <string.h>
  36. #include <math.h>
  37. #include <stdarg.h>
  38. #include <time.h>
  39. #include <ctype.h>
  40.  
  41. #if defined(MSDOS) || defined(DOS386)         /* non-blocking IO stuff */
  42. #include <io.h>
  43. #include <conio.h>
  44. #include <dos.h>
  45. #else
  46. #if defined(AMIGA_SC_6_1)
  47. #include <fcntl.h>
  48. #include <ios1.h>
  49. #include <exec/types.h>
  50. #include <exec/tasks.h>
  51. #include <dos/dos.h>
  52. #include <proto/exec.h>
  53. #include <proto/dos.h>
  54. #else
  55. #include <fcntl.h>
  56. #include <termio.h>
  57. #endif
  58. #endif
  59.  
  60. #include <setjmp.h>
  61. #define STANDARD    stderr    /* compatible with gnuplot philosophy */
  62.  
  63. #include "plot.h"
  64. #include "type.h"               /* own types */
  65. #include "matrix.h"
  66. #include "fit.h"
  67.  
  68. enum marq_res { OK, ERROR, BETTER, WORSE };
  69. typedef enum marq_res marq_res_t;
  70.  
  71.  
  72. #define INFINITY    1e30
  73. #define NEARLY_ZERO 1e-30
  74. #define DELTA        0.001
  75. #define MAX_DATA    2048
  76. #define MAX_PARAMS  32
  77. #define MAX_LAMBDA  1e20
  78. #define MAX_VALUES_PER_LINE    32
  79. #define MAX_VARLEN  32
  80. #define START_LAMBDA 0.01
  81. #if defined(MSDOS) || defined(OS2) || defined(DOS386)
  82. #define PLUSMINUS   "\xF1"                      /* plusminus sign */
  83. #else
  84. #define PLUSMINUS   "+/-"
  85. #endif
  86.  
  87. static double epsilon       = 1e-5;        /* convergence criteria */
  88.  
  89. static char fitfunction[80];
  90. static char fit_script[127];
  91. static char logfile[128]    = "fit.log";
  92. static char *fit_index_def  = "FIT_INDEX";
  93. static char *FIXED        = "# FIXED";
  94. static char *GNUFITLOG        = "FIT_LOG";
  95. static char *FITLIMIT        = "FIT_LIMIT";
  96. static char *FITSCRIPT        = "FIT_SCRIPT";
  97. static char *DEFAULT_CMD    = "replot";         /* if no fitscript spec. */
  98. static char *TEMPNAME        = "gfttemp.$$$";
  99.  
  100.  
  101. static FILE        *log_f = NULL;
  102.  
  103. static word        num_data,
  104.             num_params,
  105.             skip;
  106. static double        *fit_x,
  107.             *fit_y,
  108.             *err_data,
  109.             *a;
  110.  
  111. static struct udft_entry func;
  112.  
  113. typedef char fixstr[MAX_VARLEN+1];
  114. static fixstr        *par_name;
  115.  
  116. extern jmp_buf env;                /* from plot.c */
  117. extern char input_line[];            /* command.c */
  118.  
  119. #if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
  120. static    char buff = 0;    /* Unix, Linux, OS/2: buffer for getch(), kbhit() */
  121. #ifndef OS2
  122. #define getchx    getch    /* Unix or Linux */
  123. #endif
  124. #endif
  125.  
  126. /*****************************************************************
  127.              internal Prototypes
  128. *****************************************************************/
  129.  
  130. static marq_res_t marquardt (double a[], double **alpha, double *chisq,
  131.                  double *lambda, double *varianz);
  132. static boolean analyze (double a[], double **alpha, double beta[],
  133.                double *chisq, double *varianz);
  134. static void calculate (double *yfunc, double **dyda, double a[]);
  135. static void call_gnuplot (double *par, double *data);
  136. static void show_fit (int i, double chisq, double last_chisq, double *a,
  137.               double lambda, FILE *device);
  138. static boolean regress (double a[]);
  139. static int scan_num_entries (char *s, double vlist[]);
  140. static boolean is_variable (char *s);
  141. static void copy_max (char *d, char *s, int n);
  142. static double getdvar (char *varname);
  143. static void splitpath (char *s, char *p, char *f);
  144. static void subst (char *s, char from, char to);
  145. static void dbl_fprintf (FILE *f1, FILE *f2, char *s, ...);
  146. static boolean fit_interrupt ();
  147. static boolean is_empty (char *s);
  148.  
  149.  
  150. extern struct value *Gcomplex(struct value *v, double re, double im);
  151. extern struct value *Ginteger(struct value *a, int i);
  152. extern double real(struct value *val);
  153. extern void evaluate_at(struct at_type *at_ptr, struct value *val_ptr);
  154. extern int do_line ();
  155.  
  156. #if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
  157. /* non-blocking IO stuff */
  158. static int kbhit();
  159. static int getch();
  160. #endif
  161.  
  162. /*****************************************************************
  163.     in case of fatal errors
  164. *****************************************************************/
  165. void error_ex (char *s, ...)
  166. {
  167.     va_list p;
  168.     char    tmp[128],
  169.         *sp;
  170.  
  171.     strcpy (tmp, "          ");         /* start after GNUPLOT> */
  172.     va_start (p, s);
  173.     vsprintf (tmp+9, s, p);
  174.     va_end (p);
  175.     sp = strchr (tmp, '\0');
  176.     while ( *--sp == '\n' )
  177.     ;
  178.     strcpy (sp+1, "\n\n");              /* terminate with exactly 2 newlines */
  179.     fprintf (STANDARD, tmp);
  180.     if ( log_f ) {
  181.     fprintf (log_f, "BREAK: %s", tmp);
  182.     fclose (log_f);
  183.     log_f = NULL;
  184.     }
  185.     if ( func.at ) {
  186.     free (func.at);             /* release perm. action table */
  187.     func.at = (struct at_type *) NULL;
  188.     }
  189.  
  190. /* restore raw or cooked */
  191. #if defined(AMIGA_SC_6_1)
  192.     (void) rawcon (0);  /* Since we're returning to the command level, */
  193.                         /* we might as well assume non-raw mode. */
  194. #endif
  195.  
  196.     longjmp(env, TRUE);     /* bail out to command line */
  197. }
  198.  
  199.  
  200.  
  201. /*****************************************************************
  202.     fprintf to 2 streams at a time
  203. *****************************************************************/
  204. static void dbl_fprintf (FILE *f1, FILE *f2, char *s, ...)
  205. {
  206.     va_list p;
  207.  
  208.     va_start (p, s);
  209.     vfprintf (f1, s, p);
  210.     va_start (p, s);
  211.     vfprintf (f2, s, p);
  212.     va_end (p);
  213. }
  214.  
  215.  
  216. /*****************************************************************
  217.     Marquardt's nonlinear least squares fit
  218. *****************************************************************/
  219. static marq_res_t marquardt (double a[], double **alpha, double *chisq,
  220.                  double *lambda, double *varianz)
  221. {
  222.     int         j;
  223.     static double   *da,        /* delta-step of the parameter */
  224.             *temp_a,        /* temptative new params set   */
  225.             **one_da,
  226.             *beta,
  227.             *tmp_beta,
  228.             **tmp_alpha,
  229.             **covar,
  230.             old_chisq;
  231.  
  232.     /* Initialization when lambda < 0 */
  233.  
  234.     if ( *lambda < 0 ) {    /* Get first chi-square check */
  235.     one_da        = matr (num_params, 1);
  236.     temp_a        = vec (num_params);
  237.     beta        = vec (num_params);
  238.     tmp_beta    = vec (num_params);
  239.     da        = vec (num_params);
  240.     covar        = matr (num_params, num_params);
  241.     tmp_alpha   = matr (num_params, num_params);
  242.  
  243.     *lambda  = -*lambda;      /* make lambda positive */
  244.     return analyze (a, alpha, beta, chisq, varianz) ? OK : ERROR;
  245.     }
  246.     old_chisq = *chisq;
  247.  
  248.     for ( j=0 ; j<num_params ; j++ ) {
  249.     memcpy (covar[j], alpha[j], num_params*sizeof(double));
  250.     covar[j][j] *= 1 + *lambda;
  251.     one_da[j][0] = beta [j];
  252.     }
  253.  
  254.     solve (covar, num_params, one_da, 1);    /* Equation solution */
  255.  
  256.     for ( j=0 ; j<num_params ; j++ )
  257.     da[j] = one_da[j][0];            /* changes in paramss */
  258.  
  259.     /* once converged, free dynamic allocated vars */
  260.  
  261.     if ( *lambda == 0.0 ) {
  262.     free (beta);
  263.     free (tmp_beta);
  264.     free (da);
  265.     free (temp_a);
  266.     free_matr (one_da, num_params);
  267.     free_matr (tmp_alpha, num_params);
  268.     free_matr (covar, num_params);
  269.         return OK;
  270.     }
  271.  
  272.     /* check if trial did ameliorate sum of squares */
  273.  
  274.     for ( j=0 ; j<num_params ; j++ ) {
  275.     temp_a[j] = a[j] + da[j];
  276.     tmp_beta[j] = beta [j];
  277.     memcpy (tmp_alpha[j], alpha[j], num_params*sizeof(double));
  278.     }
  279.  
  280.     if ( !analyze (temp_a, tmp_alpha, tmp_beta, chisq, varianz) )
  281.     return ERROR;
  282.  
  283.     if ( *chisq < old_chisq ) {     /* Success, accept new solution */
  284.     if ( *lambda > 1e-20 )
  285.         *lambda /= 10;
  286.     old_chisq = *chisq;
  287.     for ( j=0 ; j<num_params ; j++ ) {
  288.         memcpy (alpha[j], tmp_alpha[j], num_params*sizeof(double));
  289.         beta[j] = tmp_beta[j];
  290.         a[j] = temp_a[j];
  291.     }
  292.     return BETTER;
  293.     }
  294.     else {                /* failure, increase lambda and return */
  295.     *lambda *= 10;
  296.     *chisq = old_chisq;
  297.     return WORSE;
  298.     }
  299. }
  300.  
  301.  
  302. /*****************************************************************
  303.     compute chi-square and numeric derivations
  304. *****************************************************************/
  305. static boolean analyze (double a[], double **alpha, double beta[],
  306.                double *chisq, double *varianz)
  307. {
  308.  
  309. /*
  310.  *  used by marquardt to evaluate the linearized fitting matrix alpha
  311.  *  and vector beta
  312.  */
  313.  
  314.     word    k, j, i;
  315.     double  wt, sig2i, dy, **dyda, *yfunc, *edp, *yp, *fp, **dr,
  316.         **ar, *dc, *bp, *ac;
  317.  
  318.     yfunc = vec (num_data);
  319.     dyda = matr (num_data, num_params);
  320.  
  321.     /* initialize alpha beta */
  322.     for ( j=0 ; j<num_params ; j++ ) {
  323.     for ( k=0 ; k<=j ; k++ )
  324.         alpha[j][k] = 0;
  325.     beta[j] = 0;
  326.     }
  327.  
  328.     *chisq = *varianz = 0;
  329.     calculate (yfunc, dyda, a);
  330.     edp = err_data;
  331.     yp = fit_y;
  332.     fp = yfunc;
  333.     dr = dyda;
  334.  
  335.     /* Summation loop over all data */
  336.     for ( i=0 ; i<num_data ; i++, dr++ ) {
  337.     /* The original code read: sig2i = 1/(*edp * *edp++); */
  338.     /* There are some compilers that evaluate the operation */
  339.     /* from right to left, although it is an error to do so. */
  340.     /* Hence the following modification: */
  341.     sig2i = 1/(*edp * *edp);
  342.     edp++;
  343.     dy = *yp++ - *fp++;
  344.     *varianz += dy*dy;
  345.     ar = alpha;
  346.         dc = *dr;
  347.     bp = beta;
  348.     for ( j=0 ; j<num_params ; j++ ) {
  349.         wt = *dc++ * sig2i;
  350.         ac = *ar++;
  351.             for ( k=0 ; k<=j ; k++ )
  352.         *ac++ += wt * (*dr)[k];
  353.         *bp++ += dy * wt;
  354.     }
  355.     *chisq += dy*dy*sig2i;            /* and find chi-square */
  356.     }
  357.  
  358.     *varianz /= num_data;
  359.     for ( j=1 ; j<num_params ; j++ )        /* fill in the symmetric side */
  360.     for ( k=0 ; k<=j-1 ; k++ )
  361.         alpha[k][j] = alpha [j][k];
  362.     free (yfunc);
  363.     free_matr (dyda, num_data);
  364.     return TRUE;
  365. }
  366.  
  367.  
  368. /*****************************************************************
  369.     compute function values and partial derivatives of chi-square
  370. *****************************************************************/
  371. static void calculate (double *yfunc, double **dyda, double a[])
  372. {
  373.     word    k, p;
  374.     double  tmp_a;
  375.     double  *tmp_low,
  376.         *tmp_high,
  377.         *tmp_pars;
  378.  
  379.     tmp_low = vec (num_data);          /* numeric derivations */
  380.     tmp_high = vec (num_data);
  381.     tmp_pars = vec (num_params);
  382.  
  383.     /* first function values */
  384.  
  385.     call_gnuplot (a, yfunc);
  386.  
  387.     /* then derivatives */
  388.  
  389.     for ( p=0 ; p<num_params ; p++ )
  390.     tmp_pars[p] = a[p];
  391.     for ( p=0 ; p<num_params ; p++ ) {
  392.     tmp_pars[p] = tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
  393. /*
  394.  *  the more exact method costs double execution time and is therefore
  395.  *  commented out here. Change if you like!
  396.  */
  397. /*      tmp_pars[p] *= 1-DELTA;
  398.     call_gnuplot (tmp_pars, tmp_low);  */
  399.  
  400.     tmp_pars[p] = tmp_a * (1+DELTA);
  401.     call_gnuplot (tmp_pars, tmp_high);
  402.     for ( k=0 ; k<num_data ; k++ )
  403.  
  404. /*            dyda[k][p] = (tmp_high[k] - tmp_low[k])/(2*tmp_a*DELTA); */
  405.  
  406.         dyda[k][p] = (tmp_high[k] - yfunc[k])/(tmp_a*DELTA);
  407.         tmp_pars[p] = a[p];
  408.     }
  409.  
  410.     free (tmp_low);
  411.     free (tmp_high);
  412.     free (tmp_pars);
  413. }
  414.  
  415.  
  416. /*****************************************************************
  417.     call internal gnuplot functions
  418. *****************************************************************/
  419. static void call_gnuplot (double *par, double *data)
  420. {
  421.     word    i;
  422.     struct value v;
  423.  
  424.     extern struct at_type *perm_at();
  425.     extern int scanner (char *);
  426.  
  427.     /* all in command.c */
  428.     extern char input_line[];
  429.     extern int num_tokens, c_token;
  430.     extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  431.     extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  432.     extern struct udft_entry *dummy_func;
  433.  
  434.     /* set parameters first */
  435.     for ( i=0 ; i<num_params ; i++ ) {
  436.     Gcomplex (&v, par[i], 0.0);
  437.     setvar (par_name[i], v);
  438.     }
  439.  
  440.     (void) strcpy (c_dummy_var[0], dummy_var[0]);
  441.  
  442.     if ( !func.at ) {            /* parse fitfunction, build action table */
  443.     strcpy (input_line, fitfunction);
  444.     num_tokens = scanner (input_line);    /* divide into tokens */
  445.     c_token = 0;
  446.     dummy_func = &func;
  447.     func.at = perm_at();            /* parse function */
  448.     }
  449.     for (i = 0; i < MAX_NUM_VAR; i++)
  450.     c_dummy_var[i][0] = '\0';               /* no dummy variables */
  451.  
  452.     dummy_func = &func;             /* by caution */
  453.     for (i = 0; i < num_data; i++) {
  454.     /* fit_index for multi-branch functions */
  455.     (void) Ginteger (&func.dummy_values[0], (i+1)*skip);
  456.     setvar (fit_index, func.dummy_values[0]);
  457.     /* calculate fit-function value */
  458.     (void) Gcomplex (&func.dummy_values[0], fit_x[i], 0.0);
  459.     evaluate_at (func.at,&v);
  460.     data[i] = real(&v);
  461.     }
  462.     (void) Ginteger (&func.dummy_values[0], 0);
  463.     setvar (fit_index, func.dummy_values[0]);
  464. }
  465.  
  466.  
  467. /*****************************************************************
  468.     getch that handles also function keys etc.
  469. *****************************************************************/
  470. #if defined(MSDOS) || defined(OS2) || defined(DOS386)
  471.                  /* if not MSDOS or OS/2, getchx not used*/
  472. int getchx ()
  473. {
  474.     register int c = getch();
  475.     if ( !c  ||  c == 0xE0 ) {
  476.     c <<= 8;
  477.     c |= getch ();
  478.     }
  479.     return c;
  480. }
  481. #endif
  482.  
  483.  
  484. /*****************************************************************
  485.     handle user interrupts during fit
  486. *****************************************************************/
  487. static boolean fit_interrupt ()
  488. {
  489.     int  c;
  490.  
  491.     while ( TRUE ) {
  492.     fprintf (STANDARD, "\n\n(S)top fit, (C)ontinue, (E)xecute:  ");
  493. #if defined(AMIGA_SC_6_1)
  494.     c = getch ();
  495. #else /* !AMIGA_SC_6_1 */
  496.     do {c = getchx ();} while ( kbhit() );
  497. #endif /* !AMIGA_SC_6_1 */
  498.         switch ( toupper (c) ) {
  499.     case 'S' :  fprintf (STANDARD, "Stop.");
  500.             return FALSE;
  501.     case 'C' :  fprintf (STANDARD, "Continue.");
  502.             return TRUE;
  503.     case 'E' :  {
  504.             int i;
  505.             struct value v;
  506.             char *tmp;
  507.  
  508.             tmp = *fit_script ? fit_script : DEFAULT_CMD;
  509.             fprintf (STANDARD, "executing: %s", tmp);
  510.             /* set parameters visible to gnuplot */
  511.             for ( i=0 ; i<num_params ; i++ ) {
  512.                 Gcomplex (&v, a[i], 0.0);
  513.                 setvar (par_name[i], v);
  514.             }
  515.             sprintf (input_line, tmp);
  516.             (void) do_line ();
  517.             }
  518.     }
  519.     }
  520. }
  521.  
  522.  
  523. /*****************************************************************
  524.     frame routine for the marquardt-fit
  525. *****************************************************************/
  526. static boolean regress (double a[])
  527. {
  528.     double    **covar,
  529.         **correl,
  530.         *dpar,
  531.         **alpha,
  532.         chisq,
  533.         last_chisq,
  534.         varianz,
  535.         lambda;
  536.     word    i, j;
  537.     int     c;
  538.     marq_res_t    res;
  539.     struct value val;
  540. #if defined(AMIGA_SC_6_1)
  541.     BPTR StdinHandle;
  542. #endif
  543.  
  544.     chisq   = last_chisq = INFINITY;
  545.     alpha   = matr (num_params, num_params);
  546.     lambda  = -START_LAMBDA;            /* use sign as flag */
  547.     i = 0;                    /* iteration counter  */
  548.  
  549.     /* Initialize internal variables and 1st chi-square check */
  550.     if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == ERROR )
  551.     error_ex ("FIT: error occured during fit");
  552.     res = BETTER;
  553.  
  554.     dbl_fprintf (STANDARD, log_f, "\nInitial set of free parameters:\n");
  555.     show_fit (i, chisq, chisq, a, lambda, STANDARD);
  556.     show_fit (i, chisq, chisq, a, lambda, log_f);
  557.  
  558.     /* MAIN FIT LOOP: do the regression iteration */
  559.  
  560.     do {
  561. #if defined(AMIGA_SC_6_1)
  562.     StdinHandle = chkufb (fileno (stdin))->ufbfh;
  563.     if (IsInteractive (StdinHandle) && 
  564.         WaitForChar (StdinHandle, 0) == DOSTRUE) {
  565.         fread (&c, 1, 1, stdin);
  566. #else /* !AMIGA_SC_6_1 */
  567.     if ( kbhit () ) {
  568.         do { getchx (); } while ( kbhit() );
  569. #endif /* !AMIGA_SC_6_1 */
  570.         show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
  571.         if ( !fit_interrupt () )           /* handle keys */
  572.         break;
  573.     }
  574.     if ( res == BETTER ) {
  575.             i++;
  576.             last_chisq = chisq;
  577.     }
  578.         if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == BETTER )
  579.         show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
  580.     }
  581.     while ( res != ERROR  &&  lambda < MAX_LAMBDA  &&
  582.         (res == WORSE  ||  (last_chisq-chisq)/chisq > epsilon) );
  583.  
  584.     dbl_fprintf (STANDARD, log_f, "\nAfter %d iterations the fit converged.\n", i);
  585.     dbl_fprintf (STANDARD, log_f, "final sum of squares residuals    : %g\n", chisq);
  586.     dbl_fprintf (STANDARD, log_f, "rel. change during last iteration : %g\n\n", (chisq-last_chisq)/chisq);
  587.  
  588.     if ( res == ERROR )
  589.     error_ex ("FIT: error occured during fit");
  590.  
  591.     /* compute errors in the parameters */
  592.  
  593.     for ( i=0 ; i<num_params; i++ )
  594.     for ( j=0 ; j<num_params; j++ )
  595.         alpha[i][j] /= varianz;
  596.  
  597.     /* get covariance-, Korrelations- and Kurvature-Matrix */
  598.     /* and errors in the parameters              */
  599.  
  600.     covar   = matr (num_params, num_params);
  601.     correl  = matr (num_params, num_params);
  602.     dpar = vec (num_params);
  603.  
  604.     inverse (alpha, covar, num_params);
  605.  
  606.     dbl_fprintf (STANDARD, log_f, "Final set of parameters            68.3%% confidence interval (one at a time)\n");
  607.     dbl_fprintf (STANDARD, log_f, "=======================            =========================================\n\n");
  608.     for ( i=0 ; i<num_params; i++ ) {
  609.     dpar[i] = sqrt (covar[i][i]);
  610.     dbl_fprintf (STANDARD, log_f, "%-15.15s = %-15g  %-3.3s %g\n",
  611.             par_name[i], a[i], PLUSMINUS, dpar[i]);
  612.     }
  613.  
  614.     dbl_fprintf (STANDARD, log_f, "\n\ncorrelation matrix of the fit parameters:\n\n");
  615.     dbl_fprintf (STANDARD, log_f, "               ");
  616.     for ( j=0 ; j<num_params ; j++ )
  617.     dbl_fprintf (STANDARD, log_f, "%-6.6s ", par_name[j]);
  618.     dbl_fprintf (STANDARD, log_f, "\n");
  619.  
  620.     for ( i=0 ; i<num_params; i++ ) {
  621.     dbl_fprintf (STANDARD, log_f, "%-15.15s", par_name[i]);
  622.     for ( j=0 ; j<num_params; j++ ) {
  623.         correl[i][j] = covar[i][j] / (dpar[i]*dpar[j]);
  624.         dbl_fprintf (STANDARD, log_f, "%6.3f ", correl[i][j]);
  625.         }
  626.     dbl_fprintf (STANDARD, log_f, "\n");
  627.     }
  628.  
  629.     /* restore last parameter's value (not done by calculate) */
  630.     Gcomplex (&val, a[num_params-1], 0.0);
  631.     setvar (par_name[num_params-1], val);
  632.  
  633.     /* call destruktor for allocated vars */
  634.     lambda = 0;
  635.     marquardt (a, alpha, &chisq, &lambda, &varianz);
  636.  
  637.     free_matr (covar, num_params);
  638.     free_matr (alpha, num_params);
  639.     free_matr (correl, num_params);
  640.     free (dpar);
  641.  
  642.     return TRUE;
  643. }
  644.  
  645.  
  646. /*****************************************************************
  647.     display actual state of the fit
  648. *****************************************************************/
  649. static void show_fit (int i, double chisq, double last_chisq, double *a,
  650.               double lambda, FILE *device)
  651. {
  652.     int k;
  653.  
  654.     fprintf (device, "\n\n\
  655. Iteration %d\n\
  656. chisquare : %-15g   relative deltachi2 : %g\n\
  657. deltachi2 : %-15g   limit for stopping : %g\n\
  658. lambda      : %g\n\nactual set of parameters\n\n",
  659. i, chisq, (chisq - last_chisq)/chisq, chisq - last_chisq, epsilon, lambda);
  660.     for ( k=0 ; k<num_params ; k++ )
  661.     fprintf (device, "%-15.15s = %g\n", par_name[k], a[k]);
  662. }
  663.  
  664.  
  665.  
  666. /*****************************************************************
  667.     is_empty: check for valid string entries
  668. *****************************************************************/
  669. static boolean is_empty (char *s)
  670. {
  671.     while ( *s == ' '  ||  *s == '\t'  ||  *s == '\n' )
  672.     s++;
  673.     return ( *s == '#'  ||  *s == '\0' );
  674. }
  675.  
  676.  
  677. /*****************************************************************
  678.     get next word of a multi-word string, advance pointer
  679. *****************************************************************/
  680. char *get_next_word (char **s, char *subst)
  681. {
  682.     char *tmp = *s;
  683.  
  684.     while ( *tmp==' ' || *tmp=='\t' || *tmp=='=' )
  685.     tmp++;
  686.     if ( *tmp=='\n' || *tmp=='\0' )                    /* not found */
  687.     return NULL;
  688.     if ( (*s = strpbrk (tmp, " =\t\n")) == NULL )
  689.     *s = tmp + strlen(tmp);
  690.     *subst = **s;
  691.     *(*s)++ = '\0';
  692.     return tmp;
  693. }
  694.  
  695.  
  696.  
  697. /*****************************************************************
  698.     get valid numerical entries
  699. *****************************************************************/
  700. static int scan_num_entries (char *s, double vlist[])
  701. {
  702.     int num = 0;
  703.     char c, *tmp;
  704.  
  705.     while ( (tmp = get_next_word (&s, &c)) != NULL
  706.                     &&  num <= MAX_VALUES_PER_LINE )
  707.     if ( !sscanf (tmp, "%lg", &vlist[++num]) )
  708.         error_ex ("syntax error in data file");
  709.     return num;
  710. }
  711.  
  712.  
  713. /*****************************************************************
  714.     check for variable identifiers
  715. *****************************************************************/
  716. static boolean is_variable (char *s)
  717. {
  718.     while ( *s ) {
  719.     if ( !isalnum(*s) && *s!='_' )
  720.         return FALSE;
  721.     s++;
  722.     }
  723.     return TRUE;
  724. }
  725.  
  726.  
  727. /*****************************************************************
  728.     strcpy but max n chars
  729. *****************************************************************/
  730. static void copy_max (char *d, char *s, int n)
  731. {
  732.     strncpy (d, s, n);
  733.     if ( strlen(s) >= n )
  734.     d[n-1] = '\0';
  735. }
  736.  
  737.  
  738. /*****************************************************************
  739.     Slow but portable implementation of strnicmp
  740. *****************************************************************/
  741. /* 
  742.  * please don't add further machine-specific macros, use -DSTRNICMP
  743.  * as compiler switch if necessary
  744.  */
  745.  
  746. #ifdef STRNICMP
  747. int strnicmp (char *s1, char *s2, int n)
  748. {
  749.     static char ss1[256], ss2[256];
  750.     register char *p;
  751.  
  752.     copy_max (ss1, s1, sizeof(ss1));
  753.     copy_max (ss2, s2, sizeof(ss2));
  754.     p = ss1-1;
  755.     while ( *++p )
  756.     toupper ( *p );
  757.     p = ss2-1;
  758.     while ( *++p )
  759.         toupper ( *p );
  760.     return strncmp (ss1, ss2, n);
  761. }
  762. #endif
  763.  
  764.  
  765. /*****************************************************************
  766.     first time settings
  767. *****************************************************************/
  768. void init_fit ()
  769. {
  770.     struct value val;
  771.  
  772.     /* index used for multi-branch functions */
  773.     fit_index = fit_index_def;
  774.     val.type = INTGR;
  775.     val.v.int_val = 0;
  776.     setvar (fit_index, val);
  777.     func.at = (struct at_type *) NULL;      /* need to parse 1 time */
  778. }
  779.  
  780.  
  781. /*****************************************************************
  782.         Set a GNUPLOT user-defined variable
  783. ******************************************************************/
  784. extern struct udvt_entry *first_udv;
  785.  
  786. void setvar (char *varname, struct value data)
  787. {
  788.     register struct udvt_entry *udv_ptr = first_udv,
  789.                    *last = first_udv;
  790.  
  791.     /* check if it's already in the table... */
  792.  
  793.     while (udv_ptr) {
  794.         last = udv_ptr;
  795.         if ( !strcmp (varname, udv_ptr->udv_name))
  796.         break;
  797.         udv_ptr = udv_ptr->next_udv;
  798.     }
  799.  
  800.     if ( !udv_ptr ) {                /* generate new entry */
  801.     last->next_udv = udv_ptr = (struct udvt_entry *)
  802.       alloc ((unsigned int)sizeof(struct udvt_entry), "setvar");
  803.     udv_ptr->next_udv = NULL;
  804.     }
  805.     copy_max (udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
  806.     udv_ptr->udv_value = data;
  807.     udv_ptr->udv_undef = FALSE;
  808. }
  809.  
  810.  
  811.  
  812. /*****************************************************************
  813.     Read INTGR Variable value, return 0 if undefined or wrong type
  814. *****************************************************************/
  815. int getivar (char *varname)
  816. {
  817.     register struct udvt_entry *udv_ptr = first_udv;
  818.  
  819.     while (udv_ptr) {
  820.         if ( !strcmp (varname, udv_ptr->udv_name))
  821.         return udv_ptr->udv_value.type == INTGR
  822.                ? udv_ptr->udv_value.v.int_val    /* valid */
  823.                : 0;                /* wrong type */
  824.             udv_ptr = udv_ptr->next_udv;
  825.     }
  826.     return 0;                        /* not in table */
  827. }
  828.  
  829.  
  830. /*****************************************************************
  831.     Read DOUBLE Variable value, return 0 if undefined or wrong type
  832. *****************************************************************/
  833. static double getdvar (char *varname)
  834. {
  835.     register struct udvt_entry *udv_ptr = first_udv;
  836.  
  837.     while (udv_ptr) {
  838.         if ( !strcmp (varname, udv_ptr->udv_name))
  839.         return udv_ptr->udv_value.type == CMPLX
  840.                ? udv_ptr->udv_value.v.cmplx_val.real    /* valid */
  841.                : 0;                   /* wrong type */
  842.             udv_ptr = udv_ptr->next_udv;
  843.     }
  844.     return 0;                         /* not in table */
  845. }
  846.  
  847.  
  848. /*****************************************************************
  849.     Split Identifier into path and filename
  850. *****************************************************************/
  851. static void splitpath (char *s, char *p, char *f)
  852. {
  853.     register char *tmp;
  854.     tmp = s + strlen(s) - 1;
  855.     while ( *tmp != '\\'  &&  *tmp != '/'  &&  *tmp != ':'  &&  tmp-s >= 0 )
  856.     tmp--;
  857.     strcpy (f, tmp+1);
  858.     strncpy (p, s, tmp-s+1);
  859.     p[tmp-s+1] = '\0';
  860. }
  861.  
  862.  
  863. /*****************************************************************
  864.     Character substitution
  865. *****************************************************************/
  866. static void subst (char *s, char from, char to)
  867. {
  868.     while ( *s ) {
  869.     if ( *s == from )
  870.         *s = to;
  871.     s++;
  872.     }
  873. }
  874.  
  875.  
  876. /*****************************************************************
  877.     write the actual parameters to start parameter file
  878. *****************************************************************/
  879. void update (char *pfile)
  880. {
  881.     char    fnam[256],
  882.         path[256],
  883.         tmpnam[256],
  884.         sstr[256],
  885.         pname[64],
  886.         tail[127],
  887.         *s = sstr,
  888. #if !defined(AMIGA_SC_6_1)
  889.         *tmpn,
  890. #endif
  891.         *tmp,
  892.         c;
  893.     FILE        *of,
  894.         *nf;
  895.     double    pval;
  896.  
  897.     /* split into path and filename */
  898.  
  899. #if defined(MSDOS) || defined(DOS386)
  900.     /* don't get problems during system() call */
  901.     subst (pfile, '/', '\\');
  902. #endif
  903.     splitpath (pfile, path, fnam);
  904.     if ( !(of = fopen (pfile, "r")) )
  905.     error_ex ("parameter file %s could not be read", pfile);
  906. #if defined(AMIGA_SC_6_1)
  907.     sprintf (tmpnam, "%sgf%08lx", path, FindTask (NULL));
  908. #else /* !AMIGA_SC_6_1 */
  909.     /* Under MSDOS rename will not work properly if TMP points to another
  910.      * drive/dir. In this case use 'path' directory for tempfile
  911.      */
  912. #if defined (MSDOS)  ||  defined (OS2)
  913.     sprintf (tmpnam, "%s%s", path, TEMPNAME);
  914. #else
  915.     tmpn = tempnam (path, "gnuft");
  916.     if ( tmpn == NULL )
  917.         error_ex ("no more unique names possible during updating");
  918.     strcpy (tmpnam, tmpn);
  919.     free (tmpn);
  920. #endif
  921. #endif /* !AMIGA_SC_6_1 */
  922.     if ( !(nf = fopen (tmpnam, "w")) )
  923.     error_ex ("could not open temporary file");
  924.     while ( TRUE ) {
  925.     if ( !fgets (s = sstr, sizeof(sstr), of) )    /* EOF found */
  926.         break;
  927.     if ( is_empty(s) ) {
  928.         fprintf (nf, s);                /* preserve comments */
  929.             continue;
  930.         }
  931.         if ( tmp = strchr (s, '#') ) {
  932.         copy_max (tail, tmp, sizeof(tail));
  933.         *tmp = '\0';
  934.     }
  935.     else
  936.         strcpy (tail, "\n");
  937.     tmp = get_next_word (&s, &c);
  938.     if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
  939.         fclose (of);
  940.         error_ex ("syntax error in parameter file %s", fnam);
  941.     }
  942.     copy_max (pname, tmp, sizeof(pname));
  943.     /* next must be '=' */
  944.     if ( c != '=' ) {
  945.         tmp = strchr (s, '=');
  946.         if ( tmp == NULL ) {
  947.         fclose (of);
  948.         error_ex ("syntax error in parameter file %s", fnam);
  949.         }
  950.         s = tmp+1;
  951.     }
  952.     tmp = get_next_word (&s, &c);
  953.     if ( !sscanf (tmp, "%lg", &pval) ) {
  954.         fclose (of);
  955.         error_ex ("syntax error in parameter file %s", fnam);
  956.     }
  957.     if ( (tmp = get_next_word (&s, &c)) != NULL ) {
  958.         fclose (of);
  959.         error_ex ("syntax error in parameter file %s", fnam);
  960.     }
  961.  
  962.         /* now modify */
  963.  
  964.     if ( (pval = getdvar (pname)) == 0 )
  965.         pval = (double) getivar (pname);
  966.  
  967.     sprintf (sstr, "%g", pval);
  968.     if ( !strchr (sstr, '.')  &&  !strchr (sstr, 'e') )
  969.         strcat (sstr, ".0");                /* assure CMPLX-type */
  970.     fprintf (nf, "%-15.15s = %-15.15s   %s", pname, sstr, tail);
  971.     }
  972.     if ( fclose (of)  ||  fclose (nf) )
  973.     error_ex ("I/O error during update");
  974.  
  975.     if ( unlink (pfile) )     /* implies write permission !!!! */
  976.       error_ex ("parameter file %s could not be deleted", pfile);
  977.  
  978. #if defined(MSDOS) || defined(OS2)
  979.     sprintf (sstr, "ren %s %s", tmpnam, fnam);
  980.     if ( system (sstr) )
  981. #else                   /* implies write permission !!!! */
  982.     if ( rename (tmpnam, pfile) )
  983. #endif
  984.       error_ex ("unable to rename %s to %s", tmpnam, pfile);
  985. }
  986.  
  987.  
  988.  
  989. /*****************************************************************
  990.     Interface to the classic gnuplot-software
  991. *****************************************************************/
  992. void do_fit (char *fitfunc, char *datafile, int xcol, int ycol, int dycol,
  993.          char *paramsfile)
  994. {
  995.     FILE    *f;
  996.     int     i, num, num_rawdata;
  997.     static char sstr[1024];
  998.     char    c,
  999.         *s = sstr,
  1000.         *tmp;
  1001.     double    vlist[MAX_VALUES_PER_LINE],
  1002.         tmpd;
  1003.     time_t    timer;
  1004.  
  1005.  
  1006. /* Amiga: force terminal to raw mode to handle user interrupts */
  1007. #if defined(AMIGA_SC_6_1)
  1008.     (void) rawcon (1);
  1009. #endif
  1010.  
  1011.     strcpy (fitfunction, fitfunc);
  1012.     if ( !(f = fopen (datafile, "r")) )
  1013.     error_ex ("datafile %s could not be read", datafile);
  1014.  
  1015.     tmpd = getdvar (FITLIMIT);          /* get epsilon if given explicitly */
  1016.     if ( tmpd < 1.0  &&  tmpd > 0.0 )
  1017.         epsilon = tmpd;
  1018.  
  1019.     *fit_script = '\0';                   /* FIT_SKIP 1 means each 2nd value */
  1020.     if ( tmp = getenv (FITSCRIPT) )
  1021.     strcpy (fit_script, tmp);
  1022.     skip = 1 + (word) getivar (FIT_SKIP);
  1023.  
  1024.     tmp = getenv (GNUFITLOG);          /* open logfile */
  1025.     if ( tmp != NULL ) {
  1026.         char *tmp2 = &tmp[strlen(tmp)-1];
  1027.         if ( *tmp2 == '/'  ||  *tmp2 == '\\' )
  1028.             sprintf (logfile, "%s%s", tmp, logfile);
  1029.         else
  1030.             strcpy (logfile, tmp);
  1031.     }
  1032.     if ( !(log_f = fopen (logfile, "a")) )
  1033.         error_ex ("could not open log-file %s", logfile);
  1034.     fprintf (log_f, "\n\n*******************************************************************************\n");
  1035.     time (&timer);
  1036.     fprintf (log_f, "%s\n\n", ctime (&timer));
  1037.     fprintf (log_f, "FIT:    data from %s\n        x = column #%d\n        y = column #%d\n",
  1038.                 datafile, xcol, ycol);
  1039.  
  1040.     fit_x = vec (MAX_DATA);         /* start with max. value */
  1041.     fit_y = vec (MAX_DATA);
  1042.  
  1043.     /* first read in experimental data */
  1044.  
  1045.     err_data = vec (MAX_DATA);
  1046.     num_rawdata = num_data = 0;
  1047.     while ( TRUE ) {
  1048.     if ( ! fgets (s = sstr, sizeof(sstr), f) )
  1049.         break;                    /* EOF found */
  1050.     if ( is_empty(s) )
  1051.         continue;
  1052.     if ( num_data==MAX_DATA ) {
  1053.         fclose (f);
  1054.         error_ex ("max. # of datapoints %d exceeded", MAX_DATA);
  1055.     }
  1056.     /* expecting n value entries (x, y (, yerror) ) */
  1057.     num = scan_num_entries (s, vlist);
  1058.     if ( xcol>num || ycol>num || dycol>num ) {
  1059.         fclose (f);
  1060.         error_ex ("columns in %s do not match specification", datafile);
  1061.     }
  1062.     if ( ++num_rawdata % skip )            /* ignore this one */
  1063.         continue;
  1064.     fit_x[num_data] = vlist[xcol];
  1065.     fit_y[num_data]   = vlist[ycol];
  1066.     err_data[num_data++] = dycol ? vlist[dycol] : 1;
  1067.     }
  1068.     fclose (f);
  1069.     /* now resize fields */
  1070.     redim_vec (&fit_x, num_data);
  1071.     redim_vec (&fit_y, num_data);
  1072.     redim_vec (&err_data, num_data);
  1073.     
  1074.     fprintf (log_f, "        #datapoints = %d", num_rawdata);
  1075.     if ( skip > 1 )
  1076.     fprintf (log_f, "   use each %d. value\n", skip);
  1077.     else
  1078.     fprintf (log_f, "\n");
  1079.     if ( dycol )
  1080.     fprintf (log_f, "        y-error = column #%d\n\n", dycol);
  1081.     else
  1082.     fprintf (log_f, "        y-errors assumed equally distributed\n\n");
  1083.     fprintf (log_f, "function used for fitting: %s\n", fitfunction);
  1084.  
  1085.     /* read in parameters */
  1086.  
  1087.     a = vec (MAX_PARAMS);
  1088.     if ((par_name = (fixstr *) malloc ((MAX_PARAMS+1)*sizeof(fixstr))) == NULL)
  1089.     error_ex ("Memory running low during fit");
  1090.     num_params = 0;
  1091.  
  1092.     if ( !strnicmp (paramsfile, "via ", 4) ) {   /* no values given */
  1093.     char *tmp2;
  1094.  
  1095.     fprintf (log_f, "use actual parameter values\n\n");
  1096.     paramsfile += 4;
  1097.     do {
  1098.         tmp = strchr (paramsfile, ',');
  1099.         if ( tmp != NULL )
  1100.         *tmp = '\0';
  1101.         tmp2 = get_next_word (¶msfile, &c);
  1102.         if ( tmp2 == NULL )
  1103.         error_ex ("no parameter specified");
  1104.         copy_max (par_name[num_params], tmp2, sizeof(fixstr));
  1105.         a[num_params] = getdvar (par_name[num_params]);
  1106.         num_params++;
  1107.         paramsfile = tmp+1;
  1108.     }
  1109.     while ( tmp );
  1110.     }
  1111.     else {
  1112.     boolean fixed;
  1113.     double    tmp_par;
  1114.  
  1115.     /* get parameters and values out of file and ignore fixed ones */
  1116.  
  1117.     fprintf (log_f, "take parameters from file: %s\n\n", paramsfile);
  1118.     if ( !(f = fopen (paramsfile, "r")) )
  1119.         error_ex ("could not read parameter-file %s", paramsfile);
  1120.     while ( TRUE ) {
  1121.         if ( !fgets (s = sstr, sizeof(sstr), f) )        /* EOF found */
  1122.         break;
  1123.         if ( tmp = strstr (s, FIXED) ) {          /* ignore fixed params */
  1124.         *tmp = '\0';
  1125.         fprintf (log_f, "FIXED:  %s\n", s);
  1126.         fixed = TRUE;
  1127.         }
  1128.         else
  1129.         fixed = FALSE;
  1130.             if ( tmp = strchr (s, '#') )
  1131.         *tmp = '\0';
  1132.         if ( is_empty(s) )
  1133.         continue;
  1134.         tmp = get_next_word (&s, &c);
  1135.         if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
  1136.         fclose (f);
  1137.         error_ex ("syntax error in parameter file %s", paramsfile);
  1138.         }
  1139.         copy_max (par_name[num_params], tmp, sizeof(fixstr));
  1140.         /* next must be '=' */
  1141.         if ( c != '=' ) {
  1142.         tmp = strchr (s, '=');
  1143.         if ( tmp == NULL ) {
  1144.             fclose (f);
  1145.             error_ex ("syntax error in parameter file %s", paramsfile);
  1146.         }
  1147.         s = tmp+1;
  1148.         }
  1149.         tmp = get_next_word (&s, &c);
  1150.         if ( !sscanf (tmp, "%lg", &tmp_par) ) {
  1151.         fclose (f);
  1152.         error_ex ("syntax error in parameter file %s", paramsfile);
  1153.         }
  1154.  
  1155.         /* make fixed params visible to GNUPLOT */
  1156.         if ( fixed ) {
  1157.         struct value v;
  1158.         Gcomplex (&v, tmp_par, 0.0);
  1159.         setvar (par_name[num_params], v);   /* use parname as temp */
  1160.         }
  1161.         else
  1162.         a[num_params++] = tmp_par;
  1163.  
  1164.         if ( (tmp = get_next_word (&s, &c)) != NULL ) {
  1165.         fclose (f);
  1166.         error_ex ("syntax error in parameter file %s", paramsfile);
  1167.         }
  1168.     }
  1169.     fclose (f);
  1170.     }
  1171.     redim_vec (&a, num_params);
  1172.     par_name = (fixstr *) realloc (par_name, (num_params+1)*sizeof(fixstr));
  1173.  
  1174.     /* avoid parameters being equal to zero */
  1175.     for ( i=0 ; i<num_params ; i++ )
  1176.     if ( a[i] == 0 )
  1177.         a[i] = NEARLY_ZERO;
  1178.  
  1179.     regress (a);
  1180.  
  1181.     fclose (log_f);
  1182.     log_f = NULL;
  1183.     free (fit_x);
  1184.     free (fit_y);
  1185.     free (err_data);
  1186.     free (a);
  1187.     free (par_name);
  1188.     if ( func.at ) {
  1189.         free (func.at);                     /* release perm. action table */
  1190.         func.at = (struct at_type *) NULL;
  1191.     }
  1192.  
  1193. /* restore raw or cooked */
  1194. #if defined(AMIGA_SC_6_1)
  1195.     (void) rawcon (0);
  1196. #endif
  1197. }
  1198.  
  1199.  
  1200. /*****************************************************************
  1201.     replace for kbhit() under UNIXish systems
  1202. *****************************************************************/
  1203. #if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
  1204. int    kbhit()
  1205. {
  1206.    int r;
  1207.    struct termio tmo, tmn;
  1208.  
  1209.    if (buff)
  1210.       return 1;
  1211.  
  1212.    ioctl(fileno(stdin),TCGETA,&tmo);
  1213.    tmn = tmo;
  1214.    tmn.c_iflag = tmn.c_lflag = 0;
  1215.    tmn.c_cc[VMIN]  = 0;      /* minimum number of characters to be read */
  1216.    tmn.c_cc[VTIME] = 1;      /* time-out in tenths of a second */
  1217.    ioctl(fileno(stdin),TCSETA,&tmn);
  1218.  
  1219.    r = read(fileno(stdin),&buff,1);
  1220.    ioctl(fileno(stdin),TCSETA,&tmo);
  1221.    return (r > 0);
  1222. }
  1223.  
  1224. int    getch()
  1225. {
  1226.    char c;
  1227.    struct termio tmo, tmn;
  1228.  
  1229.    if (buff) {
  1230.       c = buff;
  1231.       buff = 0;
  1232.       return (c & 127);
  1233.    }
  1234.  
  1235.    ioctl(fileno(stdin),TCGETA,&tmo);
  1236.    tmn = tmo;
  1237.    tmn.c_iflag = tmn.c_lflag = 0;
  1238.    tmn.c_cc[VMIN]  = 1;      /* minimum number of characters to be read */
  1239.    tmn.c_cc[VTIME] = 0;      /* time-out in tenths of a second */
  1240.    ioctl(fileno(stdin),TCSETA,&tmn);
  1241.  
  1242.    while (read(fileno(stdin),&c,1) != 1)
  1243.      ;    /* do nothing, just wait */
  1244.    ioctl(fileno(stdin),TCSETA,&tmo);
  1245.    return (c & 127);
  1246. }
  1247. #endif
  1248.  
  1249.